rm(list=ls())
setwd("YOUR WORKING DIRECTORY HERE")

library(foreign)
library(dplyr)
library(ggplot2)
library(countrycode)
library(tidyr)
library(readstata13)
library(tidyselect)
library(readxl)
library(stargazer)
library(sjPlot)
library(gridExtra)
library(grid)
library(dplyr)
library(jtools)
library(broomExtra)
library(broom.mixed)
library(ggstance)
library(lmtest)
library(car)
library(coefplot)
library(ggeffects)
library(ggpubr)

df <- read.csv("apsr_religion_legitimacy.csv")

#create logged variables
df$zpoplog <- log(df$population)
df$zpcgdlog <- log(df$gdppc)

#Figure 1
plot <- df %>%
  group_by(country) %>%
  summarize(govconf=mean(govconf, na.rm=TRUE),
            lxx=mean(lxx, na.rm=TRUE),
            mxx=mean(mxsum, na.rm=TRUE))

plot <- subset(plot, !is.na(lxx))

plot1 <- ggplot(plot, aes(reorder(country, lxx), lxx))+
  geom_col(width = 0.5)+
  coord_flip()+
  labs(title="State Support for Religion", caption="RAS Data, 1990-2014") + xlab("Country") + ylab("State Support for Religion")+
  theme(text = element_text(size = 8))

#Table 1
des <- dplyr::select(df, lxx, govconf, parliamentconf, minority, relig_import, education, wellbeing, gini, mxsum, zpoplog, zpcgdlog, polity2, zdurable)
stargazer(as.data.frame(des), type="html", out="R&T_3.doc",
          covariate.labels = c("State Support for Religion", "Confidence in Government", "Confidence in Parliament", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"))

#Table 2

summary(m1 <- lmer(govconf ~ lxx + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="R&T_3.doc")

#Figure 2
p1 <- plot_model(m1, type = "pred", terms = c("lxx"), 
                 title="Confidence in Government",
                 axis.title=c("State support for religion", "Predicted Confidence"))+
  ylim(.2, 1.5)+
  theme(text = element_text(size = 12))

p2 <- plot_model(m2, type = "pred", terms = c("lxx"), 
                 title="Confidence in Parliament",
                 axis.title=c("State support for religion", "Predicted Confidence"))+
  ylim(.2, 1.5)+
  theme(text = element_text(size = 12))

p2 <- p2 + labs(caption='Models 1 & 2, Table 2, 95% CIs, calculated using sjPlot (predictors held at mean)')

grid.arrange(p1, p2, nrow=2)

#Figure 3
df$Entanglement <- df$lxentanglement
df$Relationships <- df$lxrelationships
df$Funding <- df$lxfundingreligion
df$Enforcement <- df$lxenforcereligion

summary(m1 <- lmer(govconf ~ Entanglement + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(govconf ~ Relationships + minority + relig_import + education + wellbeing + (1 | country) + gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m3 <- lmer(govconf ~ Funding + minority + relig_import + education + wellbeing + (1 | country) + gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m4 <- lmer(govconf ~ Enforcement + minority + relig_import + education + wellbeing + (1 | country) + gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

p1 <- multiplot(m1, m2, m3, m4, 
                innerCI=0,
                outerCI=1.96,
                title="Confidence in Government",
                coefficients=c("Entanglement", "Relationships", "Funding", "Enforcement"),
                legend.position="none")+
  scale_color_manual(values=c("black","black", "black", "black"))+
  theme(text = element_text(size = 12))+
  xlim(-0.15, 0.1)

summary(m5 <- lmer(parliamentconf ~ Entanglement + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m6 <- lmer(parliamentconf ~ Relationships + minority + relig_import + education + wellbeing + (1 | country)+ gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m7 <- lmer(parliamentconf ~ Funding + minority + relig_import + education + wellbeing + (1 | country)+ gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m8 <- lmer(parliamentconf ~ Enforcement + minority + relig_import + education + wellbeing + (1 | country)+ gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

p2 <- multiplot(m5, m6, m7, m8, 
                innerCI=0,
                outerCI=1.96,
                title="Confidence in Parliament",
                coefficients=c("Entanglement", "Relationships", "Funding", "Enforcement"),
                legend.position="none")+
  scale_color_manual(values=c("black","black", "black", "black"))+
  labs(caption="Models A1-A8, Tables A1 & A2, 95% CIs, plotted using coefplot")+
  theme(text = element_text(size = 12))+
  xlim(-0.15, 0.1)

grid.arrange(p1, p2, ncol=1)

#Table 3
#majority members only
maj <- subset(df, relcat==2)
summary(m1 <- lmer(govconf ~ lxx + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=maj, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=maj, REML=FALSE))

#minority members only
min <- subset(df, relcat==0)
summary(m3 <- lmer(govconf ~ lxx + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=min, REML=FALSE))
summary(m4 <- lmer(parliamentconf ~ lxx + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=min, REML=FALSE))

stargazer(m1, m2, m3, m4,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="R&T_3.doc")

#Figure 4
summary(m1 <- lmer(govconf ~ lxx*relig_import + minority + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx*relig_import + minority + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

p1 <- ggpredict(m1, c("lxx", "relig_import"))
pp1 <- plot(p1) + 
  labs(
    x="State Support for Religion",
    y="Confidence",
    title="Confidence in Government"
  ) + guides(color=guide_legend(title="Importance of Religion"))+
  ylim(0.0, 1.5)+
  theme(text = element_text(size = 12))

p2 <- ggpredict(m2, c("lxx", "relig_import"))
pp2 <- plot(p2) + 
  labs(
    x="State Support for Religion",
    y="Confidence",
    title="Confidence in Parliament",
    caption="Models A9 & A10, Table A3, 95% CIs, calculated using sjPlot (predictors held at mean)"
  ) + guides(color=guide_legend(title="Importance of Religion"))+
  ylim(0.0, 1.5)+
  theme(text = element_text(size = 12))

ggarrange(pp1, pp2, ncol=1, nrow=2, common.legend = TRUE, legend="bottom")

######################
###APPENDIX###########
######################

#TABLE A1
summary(m1 <- lmer(govconf ~ Entanglement + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(govconf ~ Relationships + minority + relig_import + education + wellbeing + (1 | country) + gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m3 <- lmer(govconf ~ Funding + minority + relig_import + education + wellbeing + (1 | country) + gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m4 <- lmer(govconf ~ Enforcement + minority + relig_import + education + wellbeing + (1 | country) + gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

stargazer(m2, m4, m3, m1,
          dep.var.caption = "",
          covariate.labels = c("Relationships","Enforcement", "Funding", "Entanglement", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="table.doc")

#TABLE A2
summary(m5 <- lmer(parliamentconf ~ Entanglement + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m6 <- lmer(parliamentconf ~ Relationships + minority + relig_import + education + wellbeing + (1 | country)+ gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m7 <- lmer(parliamentconf ~ Funding + minority + relig_import + education + wellbeing + (1 | country)+ gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m8 <- lmer(parliamentconf ~ Enforcement + minority + relig_import + education + wellbeing + (1 | country)+ gini  + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

stargazer(m6, m8, m7, m5,
          dep.var.caption = "",
          covariate.labels = c("Relationships","Enforcement", "Funding", "Entanglement", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="table.doc")

#TABLE A3
summary(m1 <- lmer(govconf ~ lxx*relig_import + minority + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx*relig_import + minority + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Importance of Religion", "Minority Member", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="table.doc")

#TABLE A4
summary(m1 <- glm(govconf ~ lxx + minority + relig_import + education + wellbeing + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + factor(country), data=df))
summary(m2 <- glm(parliamentconf ~ lxx + minority + relig_import + education + wellbeing + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + factor(country), data=df))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="R&T_3.doc")

#TABLE A5
summary(m1 <- glm(govconf ~ lxx + minority + relig_import + education + wellbeing + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + factor(country) + factor(S002VS), data=df))
summary(m2 <- glm(parliamentconf ~ lxx + minority + relig_import + education + wellbeing + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + factor(country) + factor(S002VS), data=df))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="table.doc")

#TABLE A6
summary(m1 <- lmer(govconf ~ lxx + minority + relig_import + education + wellbeing + female + age + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx + minority + relig_import + education + wellbeing + female + age + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Female", "Age", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="R&T_3.doc")

#TABLE A7
summary(m1 <- lmer(govconf ~ lxx + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + socdisgen, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + socdisgen, data=df, REML=FALSE))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability", "Societal Discrimination"), 
          type="html", out="R&T_3.doc")

#TABLE A8
summary(m1 <- lmer(govconf ~ lxx + minority + service_freq + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx + minority + service_freq + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable, data=df, REML=FALSE))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability"), 
          type="html", out="R&T_3.doc")

#TABLE A9
df$region <- as.factor(df$region)
summary(m1 <- lmer(govconf ~ lxx + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + region, data=df, REML=FALSE))
summary(m2 <- lmer(parliamentconf ~ lxx + minority + relig_import + education + wellbeing + (1 | country) + gini + mxsum + zpoplog + zpcgdlog + polity2 + zdurable + region, data=df, REML=FALSE))

stargazer(m1, m2,
          dep.var.caption = "",
          column.labels = c("Government", "Parliament"),
          covariate.labels = c("State Support for Religion", "Minority Member", "Importance of Religion", "Education", "Wellbeing", "Gini", "Minority Discrimination", "Population (logged)", "GDP per capita (logged)", "Polity Score", "Durability", "Former USSR", "Latin America", "Sub-Saharan Africa", "Western Democracies"), 
          type="html", out="R&T_3.doc")

#TABLE A10
table <- df %>%
  group_by(country, S002)%>%
  count()
#transport this information manually to table in microsoft word

#Table A11
#manual entry of information into microsoft word
